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Ground-based atmospheric Cerenkov telescopes are proven to be effective instruments for observations of very 
high energy (VHE) 7-radiation from celestial objects. For effective use of such technique one needs detailed Monte 
Carlo simulations of 7-ray- and proton/nuclei-induced air showers in Earth atmosphere. Here we discuss in detail 
the algorithms used in the numerical code ALTAI developed particularly for the simulations of Cerenkov light 
emission from air showers of energy below 50 TeV. The specific scheme of sampling the charged particle transport 
in the atmosphere allows the performance of very fast and accurate simulations used for interpretation of the 
VHE 7-ray observations. 



1. Introduction 

Recent exciting detections and observations of 
TeV 7-ray emission from a number of galactic 
and extragalactic objects (Ong 1998, Catanese & 
Weekes 1998) have shown the high performance 
of currently operating imaging air Cerenkov tele- 
scopes (lACTs) . Several projects for future detec- 
tors have been proposed lately. The major forth- 
coming instruments, such as CANGAROO IV, 
H.E.S.S., MAGIC and VERITAS, will have sig- 
nificantly better sensitivity to 7-ray fluxes in the 
energy range from 50 GeV up to 50 TeV. 

For a lack of a coUimated test beam of TeV 
7-ray photons, the ground-based Cerenkov detec- 
tors heavily rely on the Monte Carlo simulations 
of the Cerenkov light emission from air showers 
which are used to understand the performance of 
detector. Basic parameters of the instrument (de- 
tection area, angular and energy resolution, effi- 
ciency of cosmic ray rejection, etc) can be derived 
from the simulations. The crucial point is a pre- 
cise measurement of the primary shower energy. 
For that one should include properly into the sim- 
ulations all processes of Cerenkov light emission 
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in the air shower as well as photon propagation on 
the way from the emitting shower particle to the 
telescope camera. Measurements of 7-ray fluxes, 
energy spectra, upper limits strongly depend on 
the absolute energy calibration of a telescope. 

Here we give a description of ALTAI ^ com- 
putational code developed for detailed Monte 
Carlo simulations of the Cerenkov light emission 
in TeV air showers. Among the other existing 
codes intended for such simulations, MOCCA 
(HiUas 1979), KASCADE (Kertzman & Sem- 
broski 1994), CORSIKA (Heck et al. 1997), this 
code has a distinctive advantage - it's rather high 
speed of shower simulations due to a specific algo- 
rithm used for processing the multiple scattering 
of charged shower particles. This approach does 
not consume a lot of CPU time and allows to 
perform fast and accurate simulations. Together 
with the additional routine recently developed for 
the simulations of telescope response (Hemberger 
1998) the ALTAI code was extensively used for 
production of a standard Monte Carlo database 
used in the HEGRA (High Energy Gamma Ray 
Astronomy) VHE 7-ray experiment (Konopelko 
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Figure 1. Transformation of phase coordinates 
of a charged particle at the multiple scattering 
segment. 



et al. 1999a). 

The ALTAI code consists of two major pro- 
grams which simulate the development of the 
electromagnetic (EMCCS) (hereafter we put in 
brackets the name of the corresponding subrou- 
tine of the code) and proton-nuclei (STEPAD, 
MULTIP, XPI) cascade in Earth atmosphere. 
We discuss in detail the procedure of simulating 
the electromagnetic and hadron-nuclei cascade in 
Sections 2 and 3, respectively. Section 2 also de- 
scribes the algorithms of Cerenkov light emission 
by the shower charged particle. Section 4 deals 
with nucleus- nucleus interactions. In Section 5 
we review the results of the Cerenkov light sim- 
ulations using ALTAI code compared with other 
relevant simulations and observational data. 

2. Electromagnetic cascade 

Here we overview the part of the code in- 
tended for the generation of an electromagnetic 
cascade in the atmosphere (EMCCS). Note that 
at present all cross-sections for electron and pho- 
ton interactions arc established and are well de- 
scribed in detail in the relevant energy range. 
An electromagnetic part of the code accounts 
for the following interaction processes: electron- 
positron pair production, Compton scattering for 
the primary photons and bremsstrahlung, ioniza- 
tion losses, Coulomb scattering for the primary 



electrons and positrons. 

The cross-section of e~ , e"*" pair production by 
primary photon was calculated according to the 
Bethe-Heitler formula taken from Motz et al. 
(1964) (MPRH). The emission angle of e-,e+ 
components, defined with respect to the direc- 
tion of initial photon, was sampled using the 
Bethe distribution (Motz et al. 1964). We used 
the Klcin-Nishina formula (see e.g., Gaisser 1990) 
for the cross-section of Compton photon scatter- 
ing (COMPH). The dependence of photon cross- 
sections on energy and the atomic number of a 
medium was derived from data taken from Storm 
& Israel (1973) and Hubbel (1969) (FMPM). 

We simulated the bremsstrahlung interaction 
process (MRDH) according to the Bethe-Heitler 
formula (see Koch & Motz, 1959). The angle of 
emitted photon was sampled from the Schicff dis- 
tribution (see Koch & Motz, 1959). The Bethe- 
Bloch formula, with the corrections for the den- 
sity effect (Stcrnhcimcr 1952), was used for sim- 
ulation of the mean ionization losses (MIONH) 
whereas for the differential cross-section of e~ , e+ 
ionization collisions we used the MoUer formula 
(MoUer 1932). Note that in the electromagnetic 
shower positrons were treated similarly to elec- 
trons. 

2.1. Multiple scattering 

The rclativistic electron suffers an cnormoiis 
number of interactions along its path length in the 
matter. Apparently, a direct simulation of all in- 
teractions is not possible. That is why all numer- 
ical codes for the shower simulation use a specific 
technique for grouping electron interactions (see 
for review Berger 1963, Akkerman 1991, Gaisser 
1990). In this approach one can simulate only the 
so-called catastrophic e~ , e+ interactions. These 
interactions lead to emission of photon/electron 
of a relatively high energy, E > Q, which exceeds 
the intermediate energy threshold of Q. The same 
approach was used for the simulation of the e~ , e"*" 
ionization process (Q = Qi) and bremsstrahlung 
(Q = Qr). The emitted electron must have a ki- 
netic energy above a certain intermediate energy 
threshold, T > Qj. 

A random path length between neighboring 
catastrophic collisions was divided into small seg- 
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ments of Al < Sgr/cm^. We sampled the phase 
coordinates of a charged particle at the end of 
the segment Al according to the cumulative ef- 
fect of all low energy collisions along the segment 
(SGMEH, SGMQH). Thus we simulated at the 
end of each segment the loss of the e^,e+ ki- 
netic energy, AT = Ti — T2, the scattering angle, 
6 = arccos(Oif22), the azimuth angle of scattered 
particle, Lp, the longitudinal Z and lateral p dis- 
placement of initial particle while passing over the 
segment (see Figure 1). Appropriate transport 
equations were used in order to derive the prob- 
ability distributions for the phase coordinates at 
the end of a segment. We have obtained the an- 
alytical solutions for these distributions. 

2.1.1. Energy losses 

The energy losses of a charged particle, AE, in 
the segment, Al, were sampled according to the 
Landau- Vavilov formula (Landau 1944, Vavilov 
1957). Plyasheshnikov & Kolchuzkin (1975) have 
tabulated this formula for a specific conditions 
of shower simulations. The distribution of the 
energy losses at the end of a multiple-scattering 
segment was described as follows: 
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is the mean energy loss per unit path length 
due to the inelastic collisions of a charged par- 
ticle by production of low energy secondaries. 



Wi(Q),Wr(Q) are the differential cross-sections 
for the ionization (I) and bremsstrahlung (R) in- 
teractions. 

The standard inverse function technique was 
applied in order to simulate the energy loss 
AT according to the two-dimensional distribu- 
tion /(yU, A) (Eqn. (1)). For that one needs to 
compute the tables of the function A = A(/i, a), 
which is inverse to the integral distribution 
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These tables can be found, e.g., in Akimov et.al. 
(1981) where A = A(/x, a) was tabulated over a 
two-dimensional lattice with the steps of 0.025 
and 0.01 over a and /x, respectively. First, it is 
necessary to generate the random number a uni- 
formly distributed in the interval (0,1) and cal- 
culate parameter ij, according to Eqn. (1). After 
that one should interpolate A = A(yu, a) using the 
above mentioned tables and finally determine AT 



as 



AT = P1AI + XQi 
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fi^/^-X, /i < 0.1; 
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Formula (5) was derived using two asymptotics 
of initial distribution (1) for ^ = and fi = 00, 
which were sewn up at /i = 0.1. The subrou- 
tine ELSH makes calculations according to these 
algorithms. 

2.1.2. Angular deflection 

We simulated the angle of multiple scattering 
using the Moliere theory (Moliere 1948, Bethe 
1953). In addition we have improved the Moliere 
distribution by taking into account the energy 
losses at the multiple scattering segment. This 
distribution may be described as follows: 



p(e)ede = 
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To=Ti-pi{T,)Al. 



(7) 



Parameter B was determined by resolving 
the transcendental equation _B — Ini? = 
ln(0.856(xc/Xa)^)- Functions /„(n = 0, 1, 2) were 
calculated as 
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Quantities Xc and Xa are closely related to the 
differential cross-section for Coulomb scattering 
(see formula (10) for definition of the Coulomb 
cross-section). 

For a fixed segment length Al and kinetic en- 
ergy of a particle T one can determine in a few it- 
erations the parameter B, which defines the shape 
of the distribution P(©). Such calculations were 
done by use of subroutine ANGMH. 

2.1.3. Space displacement 

At the multiple scattering segments we simu- 
lated the longitudinal displacement Z of a charged 
particle. For that we used the Yang-Spenser dis- 
tribution (Yang 1951, Spencer & Coune 1962). 
This distribution was adapted for the air shower 
simulations by Plyasheshnikov & Kolchuzhkin 
(1975). This distribution may be represented by 
the following expression 
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Here Wc(T, ^) is a cross-section of the Coulomb 
scattering which determines parameters Xa and 

Xc- 

To simulate the longitudinal displacement Z we 
use approach similar to that used in Section 2.1.1 



for simulation of the energy loss AT. It is based 
on the interpolation of the two-dimensional func- 
tion ^ = ^(i^, a), which is the inverse function to 
the integral distribution 
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One can find these tables, e.g., in Akimov et.al 
(1981). This approach includes (i) sampling of 
the random number a, (ii) calculation of v on the 
basis of quantities Ti and Al, (iii) determination 
of ^ on the basis of the two-dimensional inter- 
polation using the above mentioned tables and, 
finally, calculation of Z according to the formula 



Z = Al-£,r}Af. 



(12) 



For the lateral displacement of a charged par- 
ticles at the multiple scattering segment we de- 
fined the corresponding probability distribution 
using Fermi formula (see e.g., Kolchuzhkin & 
Plyasheshnikov 1975). This distribution was de- 
rived by solving the transport equations in the 
Focker-Planck approximation where the collision 
integral corresponding to the Coulomb scatter- 
ing was replaced by a second order differen- 
tial operator (see e.g., Kolchuzhkin & Uchaikin 
1978). For the fixed angle of multiple-scattering 
= (0rr,0i/) = (0cos</?, 8sin(^) this distribu- 
tion was as follows 



-P(V') = — exp 
7r7 
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V = p/Al, 



(13) 



where the parameter 7 = rjAl defines the width of 
the distribution. As was shown by Kolchuzhkin & 
Plyasheshnikov (1975) more accurate numerical 
solution of the transport equations reduces by a 
factor of 1.5 the width of the final distribution. 
This difference may be corrected by introducing 
another definition for the parameter 7: 



7 = x; 



(14) 



where B is determined from B — InB = 
^'"■((Xc/Xa)^/1.80). The complete procedure for 
charged particle transport in multi-dimensional 
phase space of energy, angular and space coor- 
dinates was defined by a few parameters: two 
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threshold energies Qi, Qr, and the length of a 
multiple-scattering segment Al. The extensive 
test calculations using the ALTAI code revealed 
the optimum values of parameters which allow to 
perform rather fast shower simulations without 
introducing systematic errors. Thus wc used the 
following parameters for the calculational proce- 
dure: Qi ~ 0.5 To, Qr ~ 0.1 To, where To ~ 
20 MeV is the threshold energy for Cerenkov light 
emission by cascade electrons in air. 

Apparently the length of a multiple scattering 
segment, Al, is one of the basic parameters of this 
method. By use of rather small segments one can 
reduce systematic error introduced by approxima- 
tions of analytical solutions for the phase trans- 
formations at the multiple scattering segment. 
On the other hand this may slow down the proce- 
dure of shower simulation. We found that the op- 
timum length of the multiple-scattering segment 
is within the range of ~ 1 — 5 gr/cm^. Very accu- 
rate analytical solutions of a multiple-scattering 
process for a e~,e+ transport allow us to use 
in calculations such segment length without in- 
troducing systematic errors in three-dimensional 
shower development. 

In comparison to other computational codes we 
used more accurate distributions derived analyti- 
cally from the multiple scattering theory. On the 
contrary the standard approach used for exam- 
ple in EGS-IV code does not include fluctuations 
of the energy losses at the multiple-scattering seg- 
ment. In addition, in the EGS-IV code the lateral 
displacement of charged particle at the multiple- 
scattering segment is taken as p = and cor- 
respondingly the longitudinal displacement is of 
z = Al. All this necessitates a small segment 
length and makes significantly more time consum- 
ing the shower simulations. Thus for the simula- 
tions of Cerenkov light from the TcV air shower 
ALTAI code is faster by as much as a few times 
compared with the EGS-IV. 



2.2. Emission of Cerenkov Light 

A charged particle (e, /x) in an air shower can 
emit Cerenkov light when its energy exceeds a 
certain threshold energy E > Eth- The threshold 
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Figure 2. Geometry of propagation of Cerenkov 
light photons. 



energy is determined as (Prank & Tamm 1937) 
= (15) 



^2(^-1) 

where n is the refraction index in air, mo is the 
particle mass. Thus at sea level the energy thresh- 
old is Eth -20 MeV for electron and -4 GeV for 
muons. To describe the altitude dependence of 
the refraction index we used the following expres- 
sion (Beliaev et al. 1980) 

n = l + 77, r]{H)=r,o/Po-p{H) (16) 

where 770 = 2.9 ■ lO"'', po = 1.22 ■ IQ-^g/cm^, 
and p{H) is the air density at a height H above 
the sea. The model of a standard atmosphere 
(Elterman 1968) was used in the simulations. 

According to Frank & Tamm (1937) the mean 
number of Cerenkon photons emitted in a 1 cm 
pathlength by electron is described as 



dQ 
dl 



27ra( 



1 



sin^ 6*0 = 1 - 



A2' 
1 



sin^ 6 c 



(1 - {mc^/Ey)n^ 



(17) 



6 



where a is a fine structure constant a = 1/137 
and — 1° is the opening angle of the Cerenkov 
light cone. The spectral region of emission is de- 
fined by wavelengths Ai and A2. In simulations 
we included sampling of the Cerenkov light atten- 
uation in the atmosphere due to the Raleigh scat- 
tering, ozone and aerosol absorption. The cross- 
sections for these processes were calculated using 
the data of DriscoU & Vaughan (1978). 

In the shower simulation procedure (PARCHR, 
GSTCHR) (Konopelko 1990), we first define the 
total number of Cerenkov photons < > emit- 
ted at the multiple scattering segment Al. As 
was discussed above we used rather small mul- 
tiple scattering segments. Therefore, wc can as- 
sume that all Cerenkov photons are emitted from 
the center of the multiple scattering segment (see 
Figure 1). To a rather good approximation an in- 
tersection of the Cerenkov light cone with the ob- 
servation plane forms a circle of a radius rc which 
can be calculated as 



rc = {Z- Zohs) ■ 0c 



(18) 



where Z and Z^^g are the height of Cerenkov light 
emission and the height of the observation level 
above the sea, respectively. The number of pho- 
tons hitting the telescope mirror can be calcu- 
lated as 

AL 

< fiph >= - Pi^' ^obz)) < riph > (19) 

where AL is length of the circle arc between 

points B and C (see Figure 2). p(Z,Zobs) is 
a probability of the Cerenkov light attenuation 
in the atmosphere along the way of a photon 
propagation. For specific response fimctions of 
the telescope camera one may calculate corre- 
sponding average number of photoelectrons as 

< nph.c. >= C < riph >, where ^ is a photon- 
to-photoelectrons conversion efficiency. Finally 
one can use the Poisson distribution in order to 
simulate the random number of photons (pho- 
toelectrons) hitting the telescope. The photons 
reaching the telescope mirror were uniformly dis- 
tributed over the circle arc BC (see Figure 2). 

In the latest version of the ALTAI code we save 
all, or a certain fraction, of Cerenkov photons hit- 
ting the telescope. Each photon has a set of 6 



variables, v{x, y, Ox,Oy, Zo, t). Zq is a height of the 
photon emission (in g/cm^) and f is a time of 
photon arrival to the detector. Such database of 
the photon is used for detailed simulations of the 
telescope camera response (Hemberger 1998). 

3. Hadron-nuclei cascade 

For simulations of the hadron-nuclei cascade we 
used the phenomenological model of hadron inter- 
actions (Konopelko 1990) which, in the most part, 
is based on the available accelerator data. The 
energy spectra of secondary hadrons generated in 
a pA interactions were approximated by use of 
the radial scaling model (Hillas 1979). In this 
approach the energy spectra may be presented as 

xdN/dx = F^p^g){x)Hq{x,E),x ^ E/Eo (20) 

where Eo,E are the energies of the primary and 
secondary particles, respectively. Indeces p, q de- 
note the type of primary and secondary particle 
{p,n,n~^ ,n~ jTt"). The basic functions are given 
below 



1.22(1 - a;)^-^ -I- 0.98e 



-18x 



1.3(1 



0.45' ' 

,± =F^±^^o+0.57e*(^-i), 
i/^±_(^±^o) = 1 - 0.88e-i-«^ 



(21) 



For the total cross-sections of inelastic hadron 
interactions we used the data given by Shabelski 
(1986,1987) in the following form 

' 258 < El, 
_ I 273(1 + 3.29- 10-2?/-h 

^(P,n) - S 3 88 . 10-3y2) Ei<E< E2, 

258(1 + 0.07y) E > E2 

El = 100 GeV, E2 = 10^ GeV, 

y = ln{E/ 10^ GeV) (22) 
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Table 1 



Mean number of fragments created by different nuclei (Murzin 1988) . 



Primaries | Secondaries 


Z >10 


C,N,0,F 


Li, Be, B 


He 


Z>10 


0.17 


0.29 


0.26 


1.34 


C,N,0,F 




0.11 


0.24 


1.00 


Li,Be,B 






0.15 


0.51 



The cross-sections from Eqn. (22) are measured 
in mBarn and correspond to the particle interac- 
tion with air. We assume that the mean value of 
atomic number for air nuclei is of 14.4. 

We used a simplified geometrical representation 
(Murzin 1988) for the total cross-section of the 
nucleus-nucleus interactions as follows 



R = 



+ a'J' - b 



1 + A 



1/3 



(23) 



where Aq is the effective atomic mass of air {Ao = 
14.4) and b is the effective radius of nuclei overlap- 
ping zone (5 — 1.17). aN defines the cross-section 
of hadron interaction with air. 

The transverse momenta of secondary hadrons 
were calculated according to the distribution 
given by Ranft (1972) 

J{p,n,7r±)(P±) = P±- 



1 

2D 



C 



(24) 



where the corresponding parameters (B,C,D) 
were defined for the interactions with air as de- 
scribed by Ranft et al. (1972). p± is measured in 
GeV/c. 

The ALTAI code was developed for the simula- 
tions in the energy range relevant for the very 
high energy 7-ray astronomy E < 50TeV. In 
this energy range the mean number of kaons 
(K^,K°) produced in hadron- nuclei interactions 
is very small compared with the number of emit- 
ted pions (7r='=,7r°). Besides this, in this energy 
region kaons and pions exhibit similar properties 
of inelastic interactions (see, e.g., Grishin 1982). 
Thus in the shower simulations for this restricted 
energy range we can exclude the kaon produc- 
tion process by introducing appropriate correc- 
tions for the probabilities of pion production in 
the (p, n,7r^) X inelastic collisions. 



For each inelastic hadron interaction (MUL- 
TIP) we sampled first the energy of the so-called 
leading particle (p, n,7r^). This particle carries 
away the bulk of the primary energy. We as- 
sume that the energy of a leading particle E is 
uniformly distributed within the range (0, Eq) 
where Eq is the energy of a primary interact- 
ing particle. Note that the type of the leading 
particle was selected randomly amongst proton 
and neutron for interaction (j),n) + A {p,n), 
and amongst charged pion or neutral pion for in- 
teraction (7r=^,7r°) + A ^ (7r=^,7r°). Thus the 
total inelasticity coefficient was determined as 
K ~ 1 — E / Eq. At the next step we sampled 
multiple production of other secondary particles 
(tt^, n°). The type of secondary particle was ran- 
domized assuming that 7r° pions forms, on aver- 
age, one third of all pions generated in hadron in- 
teractions. The energy of secondary particle was 
simulated using the spectra given by Eqn. (po|). 
We have completed the production of secondary 
particles when their total energy exceeds the en- 
ergy of a primary particle. The energy of the last 
simulated particle was renormalized in order to 
conserve the total energy in each inelastic inter- 
action. 

The transverse momentum of the secondaries 
was simulated using Eqn. (|2^). By the renormal- 
ization of the transverse momentum of leading 
particles we allowed to fulfill the momentum con- 
servation for each hadron interaction without dis- 
tortion of stochastic properties of a multi-particle 
production process. The test calculations have 
shown that such an approach describes very well 
the initial inclusive spectra of the secondary par- 
ticles if the number of secondaries is relatively 
high (^ 1). Note that as almost identical algo- 
rithm was developed and tested by Barashenkov 
& Toneev (1972). 
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Note that we have included in the code the 
propagation of a single muon generated in the 
hadron-nuclei cascade due to the tt^ decay 
(SGG,TRMUON). The emission of the Cerenkov 
light from the muons was included in simulations 
according to the scheme described above. 

4. Nucleus-nucleus interactions 

We have implemented in the ALTAI code a 
model of independent nucleon interactions of col- 
liding nuclei with nuclei fragmentation included 
(TAFRAC, TRAN, FRAG, OVERLAP). In this 
approach we assume that all nucleons of the pro- 
jectile nucleus have the same energy determined 
as Eo/Ap. Eq and Ap are the energy and atomic 
number of projectile nucleus, respectively. All 
nucleons of the projectile nucleus which overlap 
with the target nucleus interacted independently 
with each other. The non-overlapping part of the 
projectile nucleus decayed into individual nucle- 
ons and heavier fragments. The energy of frag- 
ment with the atomic number A is defined as 
Ef = {Eq ■ A)/Ap. We simulated a random num- 
ber of a fragments according to the probabilities 
of different channels summarized in Table 1. 



a. u. 0.3 



Gamma Rays 




1.5 2 2.5 3 
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Figure 3. Distribution of a shape parame- 
ter (mean scaled Width) of Cerenkov light im- 
ages simulated with ALTAI code (histograms) 
and recorded with HEGRA system of lACTs 
(crosses). 



5. Comparison with other codes and data 

An overview of Monte Carlo results on lat- 
eral, temporal, and angular characteristics of the 
Cerenkov light in air showers of 10 GeV - 1 TeV 
was recently given by Konopelko (1997). Most of 
the Cerenkov light characteristics calculated with 
the ALTAI code reproduce rather well the results 
obtained with the MOCCA code (Hillas 1996). 
The stereoscopic observations of BL Lac object 
Mkn 501 in 1997 with the HEGRA system of 
imaging air Cerenkov telescopes allowed the first 
measurements of the parameters of the Cerenkov 
light images from the 7-ray-induced air showers. 
The HEGRA data are in excellent agreement with 
the results of Monte Carlo simulations obtained 
with ALTAI code (Aharonian et al. 1999a). Us- 
ing the Mkn 501 data sample comprising 38,000 
7-ray events we have tested in great detail the pa- 
rameters of Cerenkov light image orientation and 
shape (Konopelko et al. 1999a) (see Figure 3) as 



predicted by the simulations. 

The HEGRA stereoscopic system of 5 lACTs 
was used for the measurements of the lateral 
distribution of the Cerenkov light in the 7-ray- 
induced air showers. These measurements have 
been compared again with the Monte Carlo sim- 
ulations using the ALTAI code (Aharonian et al. 
(1998)). The simulations using the ALTAI code 
reproduce very well the measured lateral distribu- 
tions of the Cerenkov light. Note that the mea- 
sured shape of the Cerenkov light lateral distribu- 
tions is almost independent of the detector simu- 
lation procedure but strongly depends on the de- 
velopment in space of a multi-TeV 7-ray shower 
in the atmosphere. 

We made several comparisons of the shape 
and size (total number of photoelectrons) of the 
Cerenkov light images for the proton- and nuclei- 
induced air shower calculated with the ALTAI 
and CORSIKA codes (Heck et al. 1997). The 
CORSIKA code was used with the HDPM model 
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Figure 4. Distribution of image shape parame- 
ter Width for the ALTAI simulated events (his- 
togram) and cosmic ray images recorded with 
a 10 m Whipple lACT (crosses) (Konopelko 
1999b). 



of the proton-nuclei interactions. This model is 
based on supercollider data and describes rather 
well the proton-nuclei air showers of energy be- 
low 10 TeV. In addition, the HDPM model of 
CORSIKA allows to perform relatively fast cal- 
culations with respect to other more modern 
shower generators, e.g. VENUS/QGSJET. De- 
spite the different algorithms and schemes in the 
two codes the results appeared to be almost iden- 
tical (Plyasheshnikov et al. 1997). The simula- 
tions with the ALTAI and CORSIKA codes well 
reproduce each other. Note that simulations with 
the ALTAI code are essentially less time consum- 
ing. The HEGRA data for the cosmic ray air 
showers were compared with the ALTAI simula- 
tions. Good agreement between simulations and 
data provided a precise measurement of the cos- 
mic ray proton spectrum in the energy range ~1- 
5 TeV (Aharonian et al. 1999b). 

We show in Figure 3 the distribution of so- 
called mean scaled Width parameter for pri- 



mary 7-rays and cosmic rays extracted from the 
HEGRA data as well as from the Monte Carlo 
simulations. The simiilations fit very well the 
data. Recently, we compared the relevant re- 
sults of the ALTAI simulations with the cos- 
mic ray data taken with a 10 m Whipple imag- 
ing air Cerenkov telescope during the 1995/1996 
Crab Nebula observations (Konopelko 1999b). 
For these simiilations we have used the standard 
detector response functions offered by the Whip- 
ple collaboration. The simulations reproduced 
very well the shape of recorded Cerenkov light 
images (see Figure 4). 

As mentioned above the upper energy of the 
shower simulations is about 50 TeV. It is limited 
by the simplified phenomenological model of the 
hadronic cascade. Although in the case of pure 
electromagnetic showers one can extend simula- 
tions using the ALTAI code to much higher ener- 
gies without breaking any model restrictions such 
extension may need the introduction of a number 
of changes into the code for the correct treatment 
of a larger number of particles. 

6. Summary 

Here we have presented a detailed description 
of the numeric code ALTAI, which was developed 
for simulations of Cerenkov light from extensive 
air showers. This code allows to make very fast 
and accurate calculations of the response of the 
ground-based Cerenkov detectors used in VHE 7- 
ray astronomy. Although the code was designed 
for calculating the parameters and characteristics 
of the imaging air Cerenkov telescopes it can be 
rearranged with minor changes in order to simu- 
late the responses of showerfront sampling exper- 
iments like MILAGRO, or Tibet AS-7. 

We tested our computational code against the 
data taken with two currently operating detec- 
tors, the HEGRA system of imaging air Cerenkov 
telescopes and the state-of-the-art single 10 m 
Whipple telescope. These comparisons have 
proven the high precision of the simulations. The 
computational code ALTAI is an effective tool 
for producing the simulated data for VHE 7-ray 
astronomy. Note that the forthcoming imaging 
Cerenkov detectors (CANGAROO IV, H.E.S.S., 
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MAGIC, VERITAS) will need a large amount of 
simulated data. Therefore the performance of 
these instruments may benefit from the use of the 
ALTAI code. 
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